function MultiObjectiveEGO(Global)
% <algorithm> <M>
% Multi-objective efficient global optimization
% alpha --- 0.7 --- portion of samples for Kriging construction
% num_k ---   5 --- number of infill points per iteration
% H     ---  21 --- number of reference directions

%------------------------------- Reference --------------------------------
% R. Hussein, K. Deb, A Generative Kriging Surrogate Model for Constrained 
% and Unconstrained Multi-objective Optimization, in: Proc. Genet. Evol. 
% Comput. Conf. 2016, Denver, 2016: pp. 573?580. 
%--------------------------------------------------------------------------

% This function is written by Youwei He (email: 1554748356@qq.com)

    %% Parameter setting
    [alpha,num_k,H] = Global.ParameterSet(0.7,5,21);
    % parameter for AASF in equation (10)
    rho = 1e-3;
    %% Generate the initial design points
    % number of design variables
    D = Global.D;
    % number of initial design points 11*D-1
    N = Global.N;
    %% step1-1: generate initial design points using Latin Hypercube sampling
    PopDec = repmat(Global.upper-Global.lower,N,1).*lhsamp(N,D) ...
        +repmat(Global.lower,N,1);
    %% step1-2: evaluate initial design points
    Population   = INDIVIDUAL(PopDec);
    
    %% step 2: Generate the reference direction set
    [R,N_R] = UniformPoint(H,Global.M);
    R = R./sqrt(sum(R.^2,2));
    %% Optimization
    while Global.NotTermination(Population)
        %% step 2 diversity_preserver procedure: neighborhood approach
        i_direction = linspace(1,N_R,N_R);
        for i=i_direction
            for j=1:num_k
                drawnow();
                PopDecT = Population.decs;
                PopObjT = Population.objs;
                PopConT = Population.cons;% cons>=0: meet constraints
                num_con = size(PopConT,2);
                N_Pop=length(PopDecT(:,1));
                % determine whether the problem is constrained or not
                if all(all(PopConT==0)) && N_Pop==1
                    constrained=1;
                else
                    constrained=0;
                end
                %% step 3: Points_Selector procedure
                normW = sqrt(sum(R(i,:).^2,2));
                normP = sqrt(sum(PopObjT.^2,2));
                CosineP = sum((PopObjT).*R(i,:),2)./normP./normW;
                % orthogonal distance of each point to the given reference direction
                distB = normP.*sqrt(1-CosineP.^2);
                [~,indx] = sort(distB);
                N = ceil(alpha*N_Pop);
                PopObj = PopObjT(indx(1:N),:);
                PopDec = PopDecT(indx(1:N),:);
                PopCon = PopConT(indx(1:N),:);
                %% step 4-1
                index = sum(PopCon >= 0, 2) == num_con;% feasible solutions
                PopCon(index,:)=0;
                PopCon(~index,:)=-PopCon(~index,:);     
                PopObjScaled = (PopObj-repmat(min(PopObj),N,1))./(repmat(max(PopObj),N,1)-repmat(min(PopObj),N,1));  
                if constrained
                    PopConScaled = (PopCon-repmat(min(PopCon),N,1))./(repmat(max(PopCon),N,1)-repmat(min(PopCon),N,1));
                end
                %% step 4-2
                PopSmetric = zeros(N,1);
                PopSmetricT = PopObjScaled(index,:)./repmat(R(i,:),sum(index),1);
                PopSmetric(index,:)=max(PopSmetricT,[],2)+rho*sum(PopSmetricT,2);
                if constrained
                    ASF_max = max(PopSmetric(index,1));
                    PopSmetric(~index,:) = ASF_max.*ones(sum(~index),1) + sum(PopConScaled(~index,:),2);
                end
                %% step 4-3
                kriging_obj= dacefit(PopDec,PopSmetric,'regpoly0','corrgauss',1*ones(1,D),0.001*ones(1,D),1000*ones(1,D));
                f_min = min(PopSmetric);
                %% step 5: Optimization
                infill_criterion = @(x)Infill_Standard_EI(x, kriging_obj, f_min);
                best_x=rGA(infill_criterion,Global);
                % infill point too close (not used in the paper, 
                % but this happens sometimes especially for low dimensional problem)
                if min(sqrt(sum((PopDecT-best_x).^2,2)))<1E-8
                    best_x = rGA(@(x)Infill_Maximal_Distance(x, PopDecT),Global);
                end
                Population = [Population,INDIVIDUAL(best_x)];
                if size(Population,2)>=Global.evaluation
                    return
                end
            end   
        end  
    end
end

